State space applications

ARIMA(2,1,2)

The following program estimates an ARIMA(2,1,2) model with NLS and with full maximum likelihood (using the state-space object). To find the state-space representation of an unrestricted ARMA model, use the autospec menu interactively. This program calls a subroutine to get starting values for the unrestricted ARMA parameters using the fast Hannan-Rissanen algorithm. 

 
'estimation of ARIMA(2,1,2) model
'checked 3/7/2007
 
include sub_HannanRissanen.prg
 
wfcreate arima22 q 1947:1 1986:3
 
series y
y.fill 1056.5, 1063.2, 1067.1, 1080.0, 1086.8, 1106.1, 1116.3, 1125.5, 1112.4, 1105.9, 1114.3, 1103.3, 1148.2, 1181.0, 1225.3, 1260.2, 1286.6, 1320.4, 1349.8, 1356.0, 1369.2, 1365.9, 1378.2, 1406.8, 1431.4, 1444.9, 1438.2, 1426.6, 1406.8, 1401.2, 1418.0, 1438.8, 1469.6, 1485.7, 1505.5, 1518.7, 1515.7, 1522.6, 1523.7, 1540.6, 1553.3, 1552.4, 1561.5, 1537.3, 1506.1, 1514.2, 1550.0, 1586.7, 1606.4, 1637.0, 1629.5, 1643.4, 1671.6, 1666.8, 1668.4, 1654.1, 1671.3, 1692.1, 1716.3, 1754.9, 1777.9, 1796.4, 1813.1, 1810.1, 1834.6, 1860.0, 1892.5, 1906.1, 1948.7, 1965.4, 1985.2, 1993.7, 2036.9, 2066.4, 2099.3, 2147.6, 2190.1, 2195.8, 2218.3, 2229.2, 2241.8, 2255.2, 2287.7, 2300.6, 2327.3, 2366.9, 2385.3, 2383.0, 2416.5, 2419.8, 2433.2, 2423.5, 2408.6, 2406.5, 2435.8, 2413.8, 2478.6, 2478.4, 2491.1, 2491.0, 2545.6, 2595.1, 2622.1, 2671.3, 2734.0, 2741.0, 2738.3, 2762.8, 2747.4, 2755.2, 2719.3, 2695.4, 2642.7, 2669.6, 2714.9, 2752.7, 2804.4, 2816.9, 2828.6, 2856.8, 2896.0, 2942.7, 3001.8, 2994.1, 3020.5, 3115.9, 3142.6, 3181.6, 3181.7, 3178.7, 3207.4, 3201.3, 3233.4, 3157.0, 3159.1, 3199.2, 3261.1, 3250.2, 3264.6, 3219.0, 3170.4, 3179.9, 3154.5, 3159.3, 3186.6, 3258.3, 3306.4, 3365.1, 3444.7, 3487.1, 3507.4, 3520.4, 3547.0, 3567.6, 3603.8, 3622.3, 3655.9, 3661.4, 3683.3
 
'take log 1st difference
series dlogy = dlog(y)
 
'NLS/conditional MLE
equation eq1.ls(showopts,c=1e-5,m=1000) dlogy c ar(1) ar(2) ma(1) ma(2)
show eq1.output
 
'setup ARMA(2,2) in sspace
sspace ss1
ss1.append @signal dlogy = c(1) + sv1 + c(2)*sv2 + c(3)*sv3
ss1.append @state sv1 = c(5)*sv1(-1) + c(6)*sv2(-1) + [var = exp(c(4))]
ss1.append @state sv2 = sv1(-1)
ss1.append @state sv3 = sv2(-1)
 
'these starting values (same as NLS) take 64 iterations to converge
'c(1)=0.00794                    'constant
'c(2)=0.0025                     'MA(1)
'c(3)=0.0025                     'MA(2)
'c(4)=log(@var(dlogy))    'sigma
'c(5)=0.0025                     'AR(1)
'c(6)=0.0025                     'AR(2)
 
'starting values from Hannan-Rissanen
vector(6) p0
call sub_HannanRissanen(dlogy, 2, 2, 8, p0)
c(1)=p0(1)                'constant
c(2)=p0(4)                'MA(1)
c(3)=p0(5)                'MA(2)
c(4)=log(p0(6))           'sigma
c(5)=p0(2)                'AR(1)
c(6)=p0(3)                'AR(2)
 
'estimate
ss1.ml(showopts,m=1000,c=1e-5)
freeze(out1) ss1.output
show out1
 
'different starting values which fail to improve
c(1)=0                           'constant
c(2)=-1.328                      'MA(1)
c(3)=0.392                'MA(2)
c(4)=log(p0(6))           'sigma
c(5)=1.68                 'AR(1)
c(6)=-0.749                      'AR(2)
 
ss1.ml(showopts,m=1000,c=1e-5)
freeze(out2) ss1.output
show out2
 

^top

Seasonal ARIMA(0,1,1)(0,1,1) (airline model)

The following program estimates the multiplicative seasonal ARMA model, known as the "airline" model, using NLS and full MLE. This example shows how to specify a multiplicative seasonal MA model in state-space form. 

 
'estimation of airline model
'checked 3/7/2007
 
wfcreate airline m 1949:01 1960:12
 
series x
x.fill 112, 118, 132, 129, 121, 135, 148, 148, 136, 119, 104, 118, 115, 126, 141, 135, 125, 149, 170, 170, 158, 133, 114, 140, 145, 150, 178, 163, 172, 178, 199, 199, 184, 162, 146, 166, 171, 180, 193, 181, 183, 218, 230, 242, 209, 191, 172, 194, 196, 196, 236, 235, 229, 243, 264, 272, 237, 211, 180, 201, 204, 188, 235, 227, 234, 264, 302, 293, 259, 229, 203, 229, 242, 233, 267, 269, 270, 315, 364, 347, 312, 274, 237, 278, 284, 277, 317, 313, 318, 374, 413, 405, 355, 306, 271, 306, 315, 301, 356, 348, 355, 422, 465, 467, 404, 347, 305, 336, 340, 318, 362, 348, 363, 435, 491, 505, 404, 359, 310, 337, 360, 342, 406, 396, 420, 472, 548, 559, 463, 407, 362, 405, 417, 391, 419, 461, 472, 535, 622, 606, 508, 461, 390, 432
 
'airline model transformation
series y = dlog(x,1,12)
 
'NLS/conditional MLE
equation eq1.ls(showopts,c=1e-5,m=1000) y ma(1) sma(12)
show eq1.output
 
'setup ARMA(0,1)(0,1) in sspace
sspace ss1
ss1.append @signal y = sv1 + c(1)*sv2 + c(2)*sv13 + c(1)*c(2)*sv14
ss1.append @state sv1 = [var = exp(c(3))]
ss1.append @state sv2 = sv1(-1)
ss1.append @state sv3 = sv2(-1)
ss1.append @state sv4 = sv3(-1)
ss1.append @state sv5 = sv4(-1)
ss1.append @state sv6 = sv5(-1)
ss1.append @state sv7 = sv6(-1)
ss1.append @state sv8 = sv7(-1)
ss1.append @state sv9 = sv8(-1)
ss1.append @state sv10 = sv9(-1)
ss1.append @state sv11 = sv10(-1)
ss1.append @state sv12 = sv11(-1)
ss1.append @state sv13 = sv12(-1)
ss1.append @state sv14 = sv13(-1)
 
'these starting values (same as NLS) do not converge after 1000 iterations
'c(1)=0.0025                     'MA(1)
'c(2)=0.01                'MA(12)
'c(3)=log(@var(y))        'sigma
 
'use NLS estimates as starting values
c(3)=log(@var(y)) 'sigma
 
'estimate with BHHH (Marquardt does not converge after 1000 iterations)
ss1.ml(showopts,m=1000,c=1e-5,b)
show ss1.output
 

^top

Cubic spline smoothing

The following program illustrates the use of the state-space object as a filtering/smoothing routine without estimating any parameters. 

 
'cubic spline smooothing using sspace
'revised 3/7/2007
 
'change path to program path
%path = @runpath
cd %path
 
' load workfile
load nber
series y = m03054
 
smpl @first 42:06
 
'set smoothing parameter (signal-to-noise ratio)
!q = 0.005
 
'setup sspace
sspace ss1
ss1.append @signal y = sv1 + [var=1.0]
ss1.append @state sv1 = sv1(-1) + sv2(-1) + [ename=u1]
ss1.append @state sv2 = sv2(-1) + [ename=u2]
'c(1)*c(1) is the smoothing parameter (signal-to-noise ratio)
'ss1.append @evar cov(u1,u1)=c(1)*c(1)/3
'ss1.append @evar cov(u1,u2)=0.5*c(1)*c(1)
'ss1.append @evar cov(u2,u2)=c(1)*c(1)
ss1.append @evar cov(u1,u1)=!q/3
ss1.append @evar cov(u1,u2)=0.5*!q
ss1.append @evar cov(u2,u2)=!q
   
'parse (no parameters to estimate)
ss1.ml(d,m=100,c=1e-9)
 
'get smoothed estimate
ss1.makesignals(t=smooth) ys1
 
'get HP filtered series
!lambda = 14400
y.hpf(!lambda) ys2
 
'plot
graph graph1.line y ys1 ys2
graph1.setelem(1) legend(Actual)
graph1.setelem(2) legend(Cubic spline smoothing (q = !q))
graph1.setelem(3) legend(Hodrick-Prescott filter (lambda = !lambda))
graph1.legend position(.5,0) -inbox
graph1.options size(6,2)
show graph1
 

^top

QML estimation of stochastic volatility

The following program illustrates how to estimate a basic stochastic volatility model by quasi-maximum likelihood (QML). The program compares and plots the estimated log volatility series to that from a GARCH(1,1) model. Note that the standard errors reported in the state space estimates are not correct (we currently do not have an option to compute QML standard errors in statespace). 

 
'basic stochastic volatility model (12/14/2000)
'QML estimation by state space
'revised 3/7/2007
 
wfcreate u 1 1000
 
'get data
%datafile = @runpath + "svpdx.dat"
read(skiprow=1) %datafile rt
 
'estimate GARCH(1,1) model
equation eq1.arch(1,1,m=500,c=1e-5,showopts) rt c
show eq1.output
'and get conditional variance series 
smpl 946 1000
eq1.forecast rhat rhat_se gt
 
smpl @all
'transformation for QML estimation
series y = log(rt*rt)
 
'variance of log chi-square distribution
!pi = @acos(-1)
scalar s2 = 0.5*!pi*!pi
 
'specify quasi-likelihood state-space
sspace sv
sv.append @signal y = -1.27 + ht + [var=s2]
sv.append @state ht = c(1) + c(2)*ht(-1) + [var=exp(c(3))]
 
'set starting values
c(1) = 0.01
c(2) = 0.9
c(3) = 0.1
 
'estimate
show sv.ml(showopts,m=500,c=1e-7)
 
'and retrieve state series
sv.makestate(t=pred) hf
sv.makestate(t=predse) hf_se
sv.makestate(t=filt) ht
sv.makestate(t=filtse) ht_se
sv.makestate(t=smooth) hs
sv.makestate(t=smoothse) hs_se
 
'plot to compare log variance from GARCH and stochastic volatility
graph graph1.line log(gt) hf hs
graph1.options size(8,2)
graph1.addtext(0.1,-0.3) Log volatility
graph1.setelem(1) legend(GARCH(1,1))
graph1.setelem(2) legend(One-step ahead)
graph1.setelem(3) legend(Smoothed)
graph1.legend columns(1) position(0,0)
show graph1
 

^top